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?H ' ABSTRACT 

«^T Context. Star and planet formation theories predict an evolution in the density, temperature, and velocity structure as the envelope 

collapses and forms an accretion disk. While continuum emission can trace the dust evolution, spectrally resolved molecular lines are 
needed to determine the physical structure and collapse dynamics. 

Aims. The aim of this work is to model the evolution of the molecular excitation, line profiles, and related observables during low-mass 
star formation. Specifically, the signatures of disks during the deeply embedded stage (M env > M+) are investigated. 
Methods. The semi-analytic 2D axisymmetric model of Visser and collaborators has been used to describe the evolution of the density, 
stellar mass, and luminosity from the pre-stellar to the T-Tauri phase. A full radiative transfer calculation is carried out to accurately 
determine the time-dependent dust temperatures. The time- dependent CO abundance is obtained from the adsorption and thermal 
^h . desorption chemistry. Non-LTE near-IR, FIR, and submm lines of CO have been simulated at a number of time steps. 

Results. In single dish (10-20" beams), the dynamics during the collapse are best probed through highly excited 13 CO and C 18 lines, 
I which are significantly broadened by the infall process. In contrast to the dust temperature, the CO excitation temperature derived 

from submm/FIR data does not vary during the protostellar evolution, consistent with C 18 observations obtained with Herschel and 
from ground-based telescopes. The near-IR spectra provide complementary information to the submm lines by probing not only the 
cold outer envelope but also the warm inner region. The near-IR high-/ (> 8) absorption lines are particularly sensitive to the physical 
structure of the inner few AU, which does show evolution. The models indicate that observations of 13 CO and C 18 low-/ submm 
lines within a <\" (at 140 pc) beam are well suited to probe embedded disks in Stage I (M env < M*) sources, consistent with recent 
interferometric observations. High signal-to-noise ratio subarcsec resolution data with ALMA are needed to detect the presence of 
small rotationally supported disks during the Stage phase and various diagnostics are discussed. The combination of spatially and 
spectrally resolved lines with ALMA and at near-IR is a powerful method to probe the inner envelope and disk formation process 

£L^ \ during the embedded phase. 

£S| ' Key words, stars: formation - radiative transfer - accretion, accretion disks - astrochemistry - methods : numerical 

in" 

^^ ■ 1 . Introduction the excitation and kinematics of the gas on smaller scales and up 

CO ' to higher temperatures than previously possible. It is therefore 

The semi-analytical model of the collapse of protostellar en- timely to simulate the predicted molecular excitation and line 

velopes (Shu 1977; Cassen & Moosman 1981; Terebey et al. profiles within the standard picture of a collapsing envelope. The 

1984) has been used extensively to study the evolution of aim is to identify diagnostic signatures of the different physical 

/\ gas and dust from core to disk and star (Young & Evans components and stages and to provide a reference for studies of 

£h ■ 2005; Dunham et al. 2010; Visser et al. 2009). Others have ex- more complex collapse dynamics. 

- - -' plored the effects of envelope and disk parameters represen- 
tative of specific evolutionary stages on the spectral energy A problem that is very closely connected to the collapse of 
distribution (SED) and other diagnostics (Whitney et al. 2003; protostellar envelopes is the formation of accretion disks. The 
Robitaille et al. 2006, 2007; Crapsi et al. 2008; Tobin et al. presence of embedded disks was inferred from the excess of con- 
2011). These studies have focused primarily on the dust emis- tinuum emission at the smallest spatial scales through interfer- 
sion and its relation to the physical structure. On the other hand, ometric observations (eg. Keene & Masson 1990; Brown et al. 
spectroscopic observations toward young stellar objects (YSOs) 2000; J0rgensen et al. 2005, 2009). Their physical structure can 
performed by many ground-based (sub)millimeter and infrared be determined by the combined modeling of the SED and the 
telescopes also contain information on the gas structure (Evans interferometric observations (J0rgensen et al. 2005; Brinch et al. 
1999). The molecular lines are important in revealing the kine- 2007; Enoch et al. 2009). However, the excess continuum emis- 
matical information of the collapsing envelope as well as the sion at small scales can also be due to other effects of a (mag- 
physical parameters of the gas based on the molecular excita- netized) collapsing rotating envelope (pseudo-disk, Chiang et al. 
tion. The Herschel Space Observatory and the Atacama Large 2008). Therefore, resolved molecular line observations from 
Millimeter/submillimeter Array (ALMA) provide new probes of interferometers such as ALMA are needed to clearly detect 
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the presence or absence of a stable rotating embedded disk 
(Brinch et al. 2007, 2008; Lommen et al. 2008; J0rgensen et al. 
2009; Tobinetal. 2012). 

A number of previous studies have modeled the line 
profiles based on the spherically symmetry inside-out col- 
lapse scenario described by Shu (1977) (eg., Zhou et al. 1993; 
Hogerheijde & Sandell 2000; Hogerheijde 2001; Lee et al. 
2004, 2005; Evans et al. 2005). Also, numerical hydrodynam- 
ical collapse models have been coupled with chemistry and 
line radiative transfer to study the molecular line evolu- 
tion in ID (Aikawaetal. 2008) and 2D (Brinch et al. 2007; 
van Weeren et al. 2009). Best-fit collapse parameters (e.g., 
sound speed and age) are obtained but depend on the temperature 
structure and abundance profile of the model. Visser et al. (2009) 
and Visser & Dullemond (2010) developed 2D semi-analytical 
models that describe the density and velocity structure as mat- 
ter moves onto and through the forming disk. This model has 
been coupled with chemistry (Visser et al. 2009, 2011), but no 
line profiles have yet been simulated. The current paper presents 
the first study of the CO molecular line evolution within 2D disk 
formation models. CO and its isotopologs are chosen because it 
is a chemically stable molecule and readily observed. 

Observationally, most early studies of low-mass embedded 
YSOs focused on the low-/ (J u < 6) (sub-)millimeter CO 
lines in 20-3 0" beams, thus probing scales of a few thou- 
sand AU in the nearest star- forming regions (e.g., Belloche et al. 
2002; J0rgensen et al. 2002; Lee et al. 2004; Young et al. 2004; 
Crapsi et al. 2005). More recently, ground-based high-frequency 
observations of large samples up to J u =7 are becoming routinely 
available (Hogerheijde et al. 1998; van Kempen et al. 2009a,b,c) 
and Herschel-HIFI (de Graauw et al. 2010) has opened up spec- 
trally resolved observations of CO and its isotopologs up to 
J u = 16 (E u = 660 K) (Yildizetal. 2010, 2012). In addition, 
the PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) 
instruments provide spectrally unresolved CO data from J u = 4 
to J u = 50 (E u = 55 - 7300 K), revealing multiple temperature 
components (e.g., van Kempen et al. 2010; Herczeg et al. 2012; 
Goicoechea et al. 2012; Manoj et al. 2013). Although the inter- 
pretation of the higher lines requires additional physical pro- 
cesses than those considered here (Visser et al. 2012), our mod- 
els provide a reference frame within which to analyze the lower- 
/ lines. 

Complementary information on the CO excitation is ob- 
tained from near-IR (NIR) observations of the 4.7 yum fundamen- 
tal CO (v=l-0) band seen in absorption toward YSOs. Because 
the absorption probes a pencil beam line of sight toward the cen- 
tral source, the lines are more sensitive to the inner dense part 
of the envelope than the submillimeter emission line data, which 
are dominated by the outer envelope (see Fig. 1). The NIR data 
generally reveal a cold (< 30 K) and a warm (> 90 K) com- 
ponent (Mitchell et al. 1990; Boogert et al. 2002; Brittain et al. 
2005; Smith et al. 2009; Herczeg et al. 2011). The cold compo- 
nent is seen toward all YSOs in all of the isotopolog absorptions. 
However, the warm temperature varies from source to source. 
The question is whether the standard picture of collapse and disk 
formation can also reproduce these multiple temperature compo- 
nents. 

A description of the physical models and methods is given 
in Section 2. The evolution of the molecular excitation and the 
resulting far-infrared (FIR) to submm lines of the pure rotational 
transitions of CO is discussed in Section 3, whereas the NIR ro- 
vibrational transitions is presented in Section 4. The implication 
of the results and whether embedded disks can be observed dur- 
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Fig. 1. Sketch of one quadrant of an embedded protostellar sys- 
tem with a disk and envelope. The FIR/millimeter/submillimeter 
emission comes from the envelope and outflow cavity walls. On 
the other hand, the NIR absorption (gray shaded region) also 
probes the warm region close to the star through a pencil beam, 
illustrating the complementarity of the techniques. The dotted 
orange lines indicate the stellar light. 

Table 1. Parameters used for the three different evolutionary 
models. 



Model 


^0 

[Hz] 


[km s- 1 ] 


[AU] 


[10 5 yr] 


M d 
[M Q ] 


^d,out 

[AU] 


1 

2 
3 


lO" 14 
lO" 14 
lO" 13 


0.26 
0.19 
0.26 


6700 

12000 

6700 


2.5 
6.3 
2.5 


0.1 

0.2 
0.4 


50 
180 

325 



ing Stage is discussed in Section 5. The results and conclusions 
are summarized in Section 6. 



2. Method 

2.1. Physical structure 

The two-dimensional axisymmetric model of Visser et al. 
(2009), Visser & Dullemond (2010) and Visser et al. (201 1) was 
used to simulate the collapse of a rotating isothermal spheri- 
cal envelope into a pre-main sequence star with a circumstellar 
disk. The model is based on the analytical collapse solutions of 
Shu (1977, hereafter S77), Cassen & Moosman (1981, hereafter 
CM81) and Terebey et al. (1984, hereafter TSC84). The forma- 
tion and evolution of the disk follow according to the as viscos- 
ity prescription, which includes conservation of angular momen- 
tum (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974). 
The dust temperature structure (T^ ust ) is a key quantity for the 
chemical evolution, so it is calculated through full 3D continuum 
radiative transfer with RADMC3D 1 , considering the protostel- 
lar luminosity as the only heating source. The gas temperature is 
set equal to the dust temperature, which has been found to be a 
good assumption for submm molecular lines (Doty et al. 2002; 
Doty et al. 2004). 

The model is modified slightly in order to be consistent with 
observational constraints. The density at r = 1000 AU (ftiooo) 



1 www.ita.uni-heidelberg.de/~dullemond/software/ 
radmc-3d 
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should be at most 10 6 cm" 3 for envelopes around low-mass 
YSOs (J0rgensen et al. 2002; Kristensen et al. 2012), but the in- 
terpolation scheme of Visser et al. (2009) violated that criterion. 
To correct this, the CM81 and TSC84 solutions are connected 
in terms of the dimension-less variable r = £lot by interpolating 
between 100r 2 and 10r 2 for ft = 10" 13 Hz and between 100r 2 
and r 2 for Qo - 10" 14 Hz, where Qo is the envelope's initial solid 
body rotation rate. The overall collapse, the structure of the disk 
and the chemical evolution are unaffected by this modification. 
The models evolve until one accretion time, t acc = Mq/M where 
M oc c\/G with c s the initial effective sound speed of the enve- 
lope, which is tabulated in Table 1 . 

Three different sets of initial conditions are used in this 
work (Table 1), which are a subset of the conditions explored 
by Visser et al. (2009) and Visser & Dullemond (2010). Each 
model begins with a 1 M© envelope. The two variables Do and c s 
affect the final disk structure and mass. Model 1 (c s = 0.26 km 
s" 1 , Qo = 10" 14 Hz) produces a disk with a final mass (t/t^c = 1) 
of 0.1 M Q and a final radius of 50 AU. Changing D to 10" 13 
Hz (Model 3) yields the most massive and largest of the three 
disks (0.4 M©, 325 AU). Starting with a lower sound speed of 
0.19 km s" 1 (Model 2) produces a disk of intermediate mass 
and size (0.2 M , 180 AU). The initial conditions also affect 
the flattened inner envelope structure. In Model 3, the flatten- 
ing of the inner envelope extends to > 300 AU by t/t acc = 0.5 
and extends up to -1000 AU by the end of the accretion phase. 
In the other two models, the extent of the flattening is similar 
to the extent of the disk, with Model 2 showing more flatten- 
ing in the inner envelope than Model 1 . The different evolution- 
ary stages are characterized by the relative masses in the differ- 
ent components (Robitaille et al. 2006): M env » M+ (Stage 0, 
t/kcc < 0.5), M env < M* but M env > M disk (Stage 1, t/t acc > 0.5) 
and M env < ^disk (Stage 2) (see Appendix B). 

Our main purpose is to investigate how the molecular lines 
evolve during different disk formation scenarios (Stage 0/1). 
Hence, the parameters were chosen to explore the formation of 
three different disk structures. Figure 2 zooms in on the final disk 
structures for the different models at the accretion time (t = t acc ) 
where the red lines outline the disk surfaces. The angular ve- 
locities within these lines are assumed to be Keplerian while re- 
gions outside these lines are described by the analytical collaps- 
ing rotating envelope solutions (Visser et al. 2009). The density 
structure is 2D axisymmetric with a 3D velocity field within the 
envelope and the disk. 



2.2. CO abundance 

The 12 CO abundance is obtained through the adsorption and 
thermal desorption chemistry as described in Section 2.7 of 
Visser et al. (2009). The main difference is that we have used 
both the forward and backward methods (Visser et al. 2011) to 
sample the trajectories through the disk and envelope. The chem- 
istry is still solved in the forward direction. 

CO completely freezes out at T d < 18 K (Visser et al. 2009). 
For the majority of the time steps used for the molecular line 
simulations, the dust temperature is well above 18 K everywhere 
except for the outer envelope beyond -4000 AU. However, due 
to the low densities within this region, CO is still predomi- 
nantly in the gas phase. Only at early time steps, early Stage 
0, a large fraction of CO is frozen out within the inner envelope. 
Constant isotope ratios of 12 C/ 13 C = 70 and 16 0/ 18 = 540 
(Wilson & Rood 1994) are used throughout the model to com- 
pute the abundances of 13 CO and C 18 0. 
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Fig. 2. Final disk structure for the three models from Table 1 at 
t = £ a cc- Top: gas density in the inner 500 AU. The solid line 
contours mark densities of 10 6 ' 7 ' 8 ' 9 cm -3 . Middle: temperature 
structure in the inner 500 AU. The temperature contours are log- 
arithmically spaced from 10 to 320 K. The 20, 50 and 100 K 
isotherms are marked by the solid line contours. Bottom: tem- 
perature structure in the inner 200 AU. In all panels, the dotted 
lines illustrate the lines of sight at inclinations of 15°, 45° and 
75° and the red line indicates the final disk surface. 



2.3. Line radiative transfer 

A fast and accurate multi-dimensional molecular excitation and 
radiative transfer code is needed to obtain observables at a 
number of time steps. We use the escape probability method 
from Bruderer et al. (2012) (based on Takahashi et al. 1983 and 
Bruderer et al. 2010). 

The most important aspect of simulating the rotational lines 
is the gridding of the physical structure. There are three compo- 
nents in the model that require proper gridding (Fig. A.2): the 
outflow cavity, the envelope (> 6000 AU) and the disk (< 300 
AU). To resolve the steep gradients between the different re- 
gions, 15 000 - 25 000 cells are used. A detailed description of 
the grid can be found in Appendix A. The line images are ren- 
dered at a number of time steps with a resolution of 0. 1 km s" 1 to 
resolve the dynamics of the collapsing envelope and the rotating 
disk. 

In addition to pure rotational lines in the submm and FIR, 
this work also explores the evolution of the NIR fundamen- 
tal v=l-0 rovibrational absorption of CO at 4.7 /mi. RADLite 
is used to render the large number of NIR ro- vibrational spec- 
tra (Pontoppidan et al. 2009). This is done by assuming that all 
molecules are in the vibrational ground state given by the non- 
LTE calculation described above. The assumption is valid con- 
sidering that the observed NIR emission lines by Herczeg et al. 
(2011) are blue-shifted by a few km s" 1 and have broader line 
widths than that expected from r gas = T dust that is being pre- 
sented here. The observed emission toward the YSOs seem to 
originate from additional physics that are not present in our mod- 
els. Thus, for the models presented here, it is valid to consider 
that the emission component from the inner region is negligible. 
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The main collisional partners are P-H2 and 0-H2 with the 
collisional rate coefficients obtained from the LAMDA database 
(Schoier et al. 2005; Yang et al. 2010). The dust opacities used 
in our model are a distribution of silicates and graphite grains 
covered by ice mantles (Crapsi et al. 2008). Finally, the gas-to- 
dust mass ratio is set to 100. 



2.4. Comparing to observations 

The molecular lines are simulated considering a source at a dis- 
tance of 140 pc. High spatial resolution in ray tracing is needed 
for the J u > 5 lines because of the small warm emitting region. 
A constant turbulent width of b = 0.8 km s" 1 is used in addition 
to the temperature and infall broadening to be consistent with the 
observed C 18 turbulent widths toward quiescent gas surround- 
ing low-mass YSOs with beams > 9" (J0rgensen et al. 2002). 

The simulated spectral cubes are convolved with Gaussian 
beams of I", 9" and 20" with the convol routine in the MIRIAD 
data reduction package (Sault et al. 1995). The telescopes that 
observe the low-/ rotational lines (/ u < 5) typically have > 15" 
beams, as does Herschel-HIFI for the higher-/ lines with J u « 
10. A 9" beam is appropriate for observations of J u > 5 per- 
formed with, e.g., the ground-based APEX telescope at 650 GHz 
and with Herschel-PACS and HIFI for J u > 14. A 1" beam simu- 
lates interferometric observations to be performed with ALMA. 
The submm lines are rendered nearly face-on (/ ~ 5°) since this 
is the simplest geometry to quantify the disk contribution. For 
studies of the evolution of the velocity field, an inclination of 45° 
is taken. The NIR lines are analyzed for 45° and 75° inclination 
to study the different excitation conditions between lines of sight 
through the inner envelope and the disk (Fig. 2). A line of sight 
of 45° probes the inner envelope and does not go through the 
disk while a line of sight of 75° grazes the top layers of the disk. 
Both a pencil beam approximation toward the center is taken as 
well as the full RADLite simulation to study the radiative trans- 
fer effects on those lines. 



2.5. Caveats 

The synthetic CO spectra were simulated without the presence 
of fore- and background material, such as the diffuse gas of the 
large-scale cloud where these YSOs are forming. The overall 
emission of the large-scale cloud can affect the observed line 
profiles within a large beam (> 15"), in particular the J u < 4 
lines of 12 CO and 13 CO and the J u < 2 lines of C 18 0. Also not 
included in the simulations are energetic components such as 
jets, shocks, UV heating and winds. These energetics strongly 
affect the 12 CO lines, especially the intensity of the J u > 6 rota- 
tional transitions (Spaans et al. 1995; van Kempen et al. 2009a; 
Visseretal. 2012). Luminosity flares associated with episodic 
accretion events can affect the CO abundance structure as dis- 
cussed in Visser & Bergin (2012). Their effect on the evolution 
of the CO line profile and excitation is beyond the scope of this 
paper. In general, a higher CO flux is expected during an accre- 
tion burst. 



3. FIR and submm CO evolution 

The FIR and submm CO lines up to the 10-9 transition (E u 
= 290-304 K) have been simulated with i=5°. The geometrical 
effects on the line intensities ( 13 CO and C 18 0) are less than 30% 
for different inclinations and the derived excitation temperatures 
differ by less than 5%, which is smaller than the rms error on the 
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Fig. 3. Normalized C 18 line profiles as functions of evolution 
for J u = 3, 6 and 9 at t/t acc = 0.13 (red), 0.50 (green) and 0.96 
(black) for i = 5° orientation for Model 1. The lines are con- 
volved to beams of 20" (top), 9" (middle) and 1" (bottom). 
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Fig. 4. C 18 3-2 (red), 6-5 (blue) and 9-8 (black) FWHM evo- 
lution within a 9" beam for the 3 different models at / = 5° 
orientation. 



derived temperatures. On the other hand, inclination strongly af- 
fects the derived moment maps and interpretation of the velocity 
field, therefore an inclination of 45° is used in Section 3.4. 



3.1. Line profiles 

The first time step used for the molecular line simulation is at 
t = 1000 years (t/t^c ~ 10" 3 ), where the central heating is not 
yet turned on. The evolution of the mass and the effective tem- 
perature is shown in the Appendix B. The effective temperature 
during this time step is 10 K and the models are still spherically 
symmetric. The TSC84 velocity profile for this time step reaches 
a maximum radial component of 10 km s" 1 at 0.1 AU. However, 
CO is frozen out in the inner envelope, so the FIR and submm 
lines (populated up to J u < 8) show neither wing emission nor 
blue asymmetry as typical signposts of an infalling envelope. 
Instead, the CO lines probe the static outer envelope at this point, 
where the low density has prevented CO from freezing out. 

The line profiles change shape once the effective source tem- 
perature increases to above a few thousand K at t > 10 4 yr (see 
Fig. 3 and Appendix B for the line profile evolution). The line 
profiles become more asymmetric (Fig. B.2) toward higher exci- 
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tation and smaller beams, because in both cases a larger contri- 
bution from warm, infalling gas in the inner envelope and opti- 
cal depth affect the emission lines. The line widths are due to a 
combination of thermal broadening, turbulent width and velocity 
structure, but examination of models without a systematic veloc- 
ity field confirm that the increase in broadening is mostly due to 
infall (see also Lee et al. 2004). Models 2 and 3 show narrower 
lines during Stage I (t/t acc > 0.5) since they are relatively more 
rotationally dominated than in Model 1. As shown in Fig. 4, the 
higher-/ (/ u > 5) lines are broadened significantly during Stage 
0. 

The line profiles in a 1" beam depend strongly on the veloc- 
ity field on small scales and show more structure than the line 
profiles in the larger beams (Figs. 3 and B.4). The lines in a 1" 
beam are broader and show multiple velocity components (in 
particular for 12 CO and 13 CO) reflecting the complex dynamics 
due to infall and rotation plus the effects of optical depth and 
high angular resolution. A combination of CO isotopolog lines 
can provide a powerful diagnostic of the velocity and density 
structure on 10 - 1000 AU scales (see also §3.4). 

In summary, the infalling gas can significantly broaden the 
optically thin gas as shown in Fig. 4. The velocity profile in the 
inner envelope affects the broadening of the high-/ (J u > 6) 
line. These lines can be a powerful diagnostic of the velocity 
and density structure in the inner 1000 AU (< 9" beams). 

3.2. CO rotational temperature 

The observed integrated intensities (K km s" 1 ) are commonly 
analyzed by constructing a Boltzmann diagram and calculating 
the associated rotational temperature. Another useful representa- 
tion is to plot the integrated flux (J F v dv in W m~ 2 ) as a function 
of upper level rotational quantum number (spectral line energy 
distribution or SLED) to determine the / level at which the peak 
of the molecular emission occurs. The conversion between inte- 
grated intensities and integrated flux is given by 

2k 



I 



F v dv 



A 3 



da 



I T v dv. 



(1) 



Examples of single-temperature fits through the modeled 
lines are shown in Fig. 5 for Model 1. Observationally, temper- 
atures are generally measured from J u = 2 to 10 (Yildiz et al. 
2012, 2013) and from 4 to 10 (Goicoechea et al. 2012). The evo- 
lution of the rotational temperatures within a 9" beam is shown 
in Fig. 6 for the three models. At very early times, t/t acc < 0.1, a 
single excitation temperature of ~ 8-12 K characterizes the CO 
Boltzmann diagrams. After the central source turns on, the exci- 
tation temperature is nearly constant with time. In a beam of 9", 
most of the flux up to J u = 10 comes from the > 100 AU region, 
which is not necessarily the > 100 K gas (Yildiz et al. 2010) 
even for face-on orientation. Therefore, within > 9" beams, the 
observed emission is not sensitive to the warming up of material 
as the system evolves. 

A value of 46 to 65 K characterizes the 12 CO distribution for 
all models. However, 12 CO is optically thick, which drives the 
temperature to higher values due to under-estimated low -J col- 
umn densities. The C 18 and 13 CO distributions are fitted with 
similar excitation temperatures between 29 and 42 K within a 
9" beam. Because Model 3 has a higher inner envelope density, 
the derived rotational temperature at the second time step is rel- 
atively high (> 40 K) due to an optical depth effect. 

The integrated fluxes can also be used to construct CO 
SLEDs. The 12 CO SLED peaks at J u = 6 - 8 with small depen- 
dence on the heating luminosity (Fig. 7). The peak of the SLED 
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Fig. 5. Examples of single-temperature fits through the 12 CO 
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does not depend on the beam size for beams > 9". The same 
applies to the 13 CO and C 18 SLEDs which peak at J u = 4 - 5. 
There is no significant observable evolution in the SLEDs once 
the star is turned on. A passively heated system is characterized 
by such SLEDs regardless of the evolutionary state as long as 
the heating luminosity is ~ 1 - 10 L (r* ~ 1000-4500 K). 

The picture within a 1" beam is different since it resolves 
the inner envelope. The best-fit single temperature component 
depends on whether or not the disk fills a significant fraction of 
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the beam. If it does (Models 2 and 3), the C 18 and 13 CO ro- 
tational temperatures are -40-50 K. There is less spread than 
within the 9" beam due to the higher densities bringing the ex- 
citation closer to LTE. Without a significant disk contribution in 
al" beam (Model 1), the 13 CO excitation temperature is closer 
to 60 K while that of C 18 is similar to the other models at ~ 40 
K. The 12 CO rotational temperature is T rot ^ 70 K in all of the 
models. Furthermore, the CO SLED within a 1" beam peaks at 
J u = 6, 8 and 10 for C 18 0, 13 CO and 12 CO, respectively. 

A comparison with the sample presented in Yildiz et al. 
(2013) indicates that the lack of T YOt evolution is consistent with 
observations, and our predicted values for T wt match the data. 
On the other hand, the observed 13 CO T rot and SLED are gener- 
ally higher than predicted, which suggests an additional heating 
component is needed to excite the 13 CO lines. A more detailed 
comparison between model prediction and observation can be 
found in Yildiz et al. (2013). 

In summary, both rotational temperatures derived from 
Boltzmann diagrams and CO SLEDs of passively heated sys- 
tems do not evolve with time once the star is turned on. Optically 
thin 13 CO and C 18 lines are characterized by single excitation 
temperatures of the order 30-40 K within a wide range of beams 
(> 9"). The CO SLED peaks at J u = 1 ± 1 and J u = 4 - 5 for 
12 CO and other isotopologs, respectively. In 1" beam, the peak 
of the SLED shifts upward by ~ 2 / levels and a warmer rota- 
tional temperature by 10 K in the presence of a disk. 

3.3. Disk contribution to fine fluxes 

After deriving a number of observables, an important question is, 
what the fraction of the flux contributed by the disk is? This con- 
tribution can be calculated from the averaged (over line profile 
and direction) escape probability of the line emission 77^, which 
is the probability that a photon escapes both the dust and line 
absorption (Takahashi et al. 1983; Bruderer et al. 2012). The es- 
cape probabilities are used to calculate the cooling rate, r coo i 5 y, 



of each computational cell and each transition with the following 
formula 



r c0 oi,v = -rhv u iA u \nx u T] uh 
An 



(2) 



where v u i is the frequency of the line, A u \ is the Einstein A coef- 
ficient, x u is the normalized population level and n is the density 
of the molecule. The disk contribution is then the ratio of the sum 
of cooling rates from the cells in the disk compared to the cells 
within a Gaussian beam. The cooling rates are weighted with a 
Gaussian beam of size 9" or 1" while the disk emitting region 
is defined as shown in Fig. 2. The flux is then given through an 
integration of line of sight, F v = Jr c0 oi,v<i^LOS, which translates 
into the same constant factor in both disk and total fluxes. 

For a single-dish 9" beam, one needs to go to higher rota- 
tional transitions with J u > 6 at later stages to obtain a > 50% 
disk contribution, if the disk can be detected at all (Fig. B.5). 
The disk is difficult to observe directly in 12 CO and 13 CO emis- 
sion since the lines quickly become optically thick unless lines 
with J u > 10 are observed. A larger disk contribution is seen in 
the optically thin C 18 lines, but here the low absolute flux may 
become prohibitive. 

For example, the expected C 18 disk fluxes for J u > 14 
within the PACS wavelength range are < 10" 20 W m -2 which 
will take > 100 hours to detect with PACS. Herschel-RIFI is 
able to spectrally resolved the C 18 J u = 10 and 9 lines but has 
a ~ 20" beam, which lowers the disk fraction by a factor of 2 
relative to a 9" beam. Peak temperatures of only 1-8 mK dur- 
ing Stage and 2-12 mK in Stage I phase are expected, which 
are readily overwhelmed by the envelope emission. The spec- 
trally resolved C 18 spectra observed with HIFI indeed do not 
show any sign of disk emission, consistent with our predictions 
(San Jose-Garcia et al. 2013; Yildiz et al. 2012, 2013). 

A more interesting result is the disk contribution to the CO 
lines within a 1" (140 AU) region as shown in Fig. 8. The sim- 
ulations suggest a significant disk contribution (> 50%) for the 
13 CO and C 18 lines within this scale, even for the / = 1-0 
transition. For Model 3, the high disk contribution is expected 
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Fig. 10. Position- velocity slice for the 13 CO and C 18 / u = 3 and 
J u = 6 transitions along the major axis of the zeroth moment 
map. The solid black lines are the Keplerian velocity structure 
calculated from the stellar mass. 



to happen before the system enters Stage 1 (M env < M+): a sig- 
nificant disk contribution is seen even at t/t acc = 0.3. However, 
Model 2 shows a significant disk contribution only in the second 
half of the accretion phase. Model 2 shows a significant jump 
in the disk contribution between t/t acc = 0.5 and 0.75 which is 
due to significant contribution from the flattened inner envelope 
within 1" at t/t^c - 0.5. Thus, the relatively more optically thin 
CO isotopolog emission is dominated by the central rotating disk 
starting from the Stage 1 phase (t/t acc > 0.5). 

3.4. Disentangling the velocity field 

How can the rotating and infalling flattened envelope be disen- 
tangled from the Keplerian motion of the disk? The evolution of 
the rotationally dominated region is consistent with that reported 
by Brinch et al. (2008) based on more detailed hydrodynamics 
simulations. All models become rotationally dominated within 
500 AU at t/t acc >0.3. 

The presence of an embedded disk can be inferred from the 
presence of elongated 13 CO and C 18 integrated intensity maps 
(moment of Fig. 9) coupled with the moment 1 map. In the 
presence of a stable disk, the zeroth moment map is elongated 
perpendicular to the outflow axis with a double peaked structure. 
This feature is not seen in Model 1 (Fig. 9 left) since the disk 
is much smaller than the 1" beam. In the case of Model 3, the 
relatively massive disk exhibits a double peaked zeroth moment 
map which is perpendicular to the outflow direction. 

A velocity gradient is seen in the moment 1 maps for both 
Models 1 and 3. With the high resolution of the modeled spec- 
tra, it is possible to differentiate between the disk and envelope 
(Hogerheijde 2001). The presence of a velocity gradient along 
the major axis of the elongation in the moment 1 map is a clear 
sign of a stable embedded disk in the case of Model 3. In addi- 
tion, from the analysis in Sect. 3.3, we can also attribute the bulk 
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of optically thin emission in Model 3 to the disk. Meanwhile, the 
infalling rotating envelope contribution can be detected through 
the fact that the moment 1 map is not perfectly aligned but 
skewed. On the other hand, a velocity gradient without the pres- 
ence of elongation in the zeroth moment map such as in Model 
1 (Fig. 9) indicates an infalling envelope. 

Position-velocity (PV) diagrams provide another way to 
study the velocity structure. Figure 10 shows PV slices follow- 
ing the direction of the major axis of the disk (PA = 90°) in 
the integrated intensity map in the 3-2 and 6-5 lines. A sim- 
ple way to analyze PV diagrams is to separate the diagram into 
four quadrants and examine which quadrants have considerable 
emission. An infall dominated PV diagram without rotation is 
symmetric about the velocity axis with all four quadrants filled 
(Ohashi et al. 1997). The presence of rotation breaks the sym- 
metry in the PV diagram and causes the peak positions to be 
off-centered. As found previously in hydrodynamical simula- 
tions by Brinch et al. (2008), infall dominates the PV diagrams 
at early times while rotation dominates the later times. These 
features have been observed toward low-mass YSOs which sug- 
gests that generally the infalling rotating envelope dominates 
the PV diagrams (eg. Sargent & Beckwith 1987; Hayashi et al. 
1993; Saito et al. 1996; Ohashi et al. 1997; Hogerheijde 2001). 

Focusing on Model 3 at t/t acc = 0.16 (near the end of Stage 
I) as shown in Fig. 10, the C 18 PV diagrams show a clearer 
pure rotation dominated structure than the 13 CO lines. They can 
thus be used to constrain the stellar mass as long as the inclina- 
tion is known. Also, the higher-/ lines provide a better view into 
the rotating structure and should give tighter constraints on the 
stellar mass. 

Another method is to plot the velocity as function of dis- 
tance of the position of the peak intensity (Sargent & Beckwith 
1987). This is done for the red-shifted and blue-shifted compo- 
nents, plotted in Fig. 11. The method is similar to the spectroas- 
trometry technique employed in the optical and NIR observa- 
tions (Takami et al. 2001; Baines et al. 2006; Pontoppidan et al. 
2008). A Keplerian disk is characterized by a v oc r -05 with large 
contribution from the high velocity gas which is optically thin. 
As noted by Sargent & Beckwith (1987), as long as the line is 
spectrally resolved (dV = 0.1km s" 1 ), the peak position corre- 
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sponds to the maximum radius of a given velocity. On the other 
hand, the emission at low velocities (dV ~ V) is relatively more 
optically thick and is dominated by the infalling rotating enve- 
lope which peaks closer to the center (no offset). Such an anal- 
ysis can be performed directly from the interferometric data and 
is a powerful tool in searching for embedded rotationally stable 
disks out of the rotating infalling envelope. 

Why is there a difference between the red- shifted and blue- 
shifted components in Fig. 11? At high velocities, this is due 
to unresolved emission hence the different velocity components 
simply peak at the same position and the difference reflects the 
uncertainty in locating the emitting region. At low velocities, the 
optically thick infalling rotating envelope affects the peak posi- 
tions. In an infalling envelope, the blue-shifted emission is rel- 
atively more optically thin than the red-shifted emission (Evans 
1999). Such difference in the optical depth causes the red- and 
blue-shifted emissions to be asymmetric as found in the moment 
1 map. 



4. NIR CO absorption lines 

A complementary probe of the molecular excitation conditions 
is provided by the NIR CO ro- vibrational absorption lines. In 
this work, we concentrate on the v = 1 - band at 4.76 /mi. This 
absorption takes place along the line of sight through the enve- 
lope and/or disk up to where the continuum is formed. The ab- 
sorption lines are therefore computed for different inclinations of 
45° and 75°. The inclinations were chosen such that they probe 
the envelope (> 15°) and the part of the disk that is not com- 
pletely optically thick such that there is enough observable NIR 
continuum (< 75°), i.e., lines of sight that graze the disk atmo- 
sphere. Since our focus is on the excitation, the absorption lines 
have been calculated without any velocity field besides a tur- 
bulent width (Doppler b) of 0.8 km s" 1 . More importantly, the 
formation of the 4.7/mi continuum is strongly dependent on the 
inclination which affect the molecular absorption lines. 



4.1. Evolution of NIR absorption spectra 

4.1 .1 . Radiative transfer and non-LTE effects 

The line center optical depth is one of the quantities derived from 
the model that can be compared to observations. This optical 
depth can either be determined by computing the line of sight 
integrated column densities and converting this to optical depth 
or by using a full RADLite calculation. For both approaches, the 
same level populations are used, i.e., the same model for the CO 
excitation is adopted as in Section 3. In the simplest method, the 
line center optical depth is obtained from 



TO = 



c 3 9u 



87r ^bv 3 g\ 



A u iM, 



(3) 



where v is the line frequency and N\ is the lower level column 
density along the line of sight. The main difference between the 
methods is that RADLite solves the radiative transfer equation 
along the line of sight and thus accounts for continuum and 
line optical depth effects and scattered continuum photons, while 
such effects are not considered in the column density approach. 
In addition, RADLite accounts for the continuum formation, 
while a ray through the center does not. Thus, the two meth- 
ods effectively compare the total mass present along a ray and 
the observed mass as probed by the line optical depth returned 
by RADLite. The optical depth is extracted from the RADLite 
spectra using the line-to-continuum ratio at the line center. 

As discussed and illustrated extensively in Appendix C.l, 
the two approaches can give very different results, especially for 
the higher / lines for which the line optical depths can differ 
by more than an order of magnitude. The main reason is illus- 
trated in Fig. C.3, which shows the region where the NIR contin- 
uum arises. For small inclinations, the continuum is essentially 
point-like, but for higher inclinations larger radii contribute sig- 
nificantly and the continuum is no longer a point source. Thus, 
for high inclinations, the high-/ lines start absorbing off-center, 
away from warm gas in the inner few AU that are included in 
the column method. The location of the continuum is typically 
at >10 AU at i = 75° which results in two to three orders of mag- 
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nitude difference in total column density. For the case of i ~ 45°, 
the result depends on the physical structure in the inner few AU 
but the absorption can miss the warm high density region close 
to the midplane of the disk where most of the J > 10 is located 
(Fig. C.2). 

Figure C.5 and C.6 also compare LTE and non-LTE spec- 
tra. Consistent with Section 3, significant differences are found 
for the higher pure rotational levels which are sub thermally ex- 
cited. However, the radiative transfer effect is generally more 
important than the non-LTE effects for viewing angles through 
the disk/inner envelope surface. 



4.1.2. NIR spectra 

The NIR absorptions are rendered with RADLite for t/t^c > 0.5 
which starts when M env ~ M* (Stage 0/1 boundary) in order 
to have a strong 4.7/mi continuum. Figure 12 show the P and 
R branches of the CO isotopologs for Model 1 and 2 which 
show the greatest difference. In general, as shown in Fig. 12, 
the higher-/ absorptions are stronger as the system evolves into 
Stage II due to the central luminosity. For Model 1, at t/t acc = 
0.5, the 12 CO absorptions can be seen up to / = 16, and up to 
J = 9 for 13 CO and 7 for C 18 0. The number of absorption lines 
that are above r y = 3 x 10" 3 (the 3<x observational limit with 
instruments such as VLT-CRIRES) doubles at the end of the ac- 
cretion phase and the highest observable / shifts from 16 to > 30 
for 12 CO and up to ~ 15-20 for the isotopologs. 

At early times, a similar number of lines are found for Model 
3 but their optical depths are lower due to the difference in the 
inner envelope structure. On the other hand, the number of lines 
in Model 2 is much higher than in the other two models due to 
the lower continuum optical depth in the inner envelope. Thus, 
the absorbing column starts deeper than in the other two models, 
which results in stronger absorption features and an increased 
number of detectable high-/ lines, up to / > 40 for 12 CO, / ~ 25 
for 13 CO and J ~ 15 for C 18 (Fig. 12). Thus, in principle the 
appearance of the NIR spectrum could be a sensitive probe of 
the inner envelope structure. In practice, the 12 CO absorption 
will be affected by outflows and winds so the 13 CO and C 18 
isotopologs are most useful for this purpose. 



4.2. Rotational temperatures 

A directly related observable is the rotational temperature de- 
rived from a Boltzmann diagram. Thus, the RADLite NIR spec- 
tra need to be converted into column densities. For this conver- 
sion, we use a standard curve-of-growth analysis (Spitzer 1989) 
with b = 0.8 km s" 1 and oscillator strengths derived from the 
Einstein A u \ coefficients. By using the curve-of-growth method, 
it is assumed that the spectral lines are unresolved. 

An example of a Boltzmann diagram is shown in Fig. C.l 
together with a two-temperature fit, as commonly done in ob- 
servations (Mitchell et al. 1990; Smith et al. 2009; Herczeg et al. 
201 1). The rotational temperature of the cold component is ob- 
tained from fitting the P(l) to P(4) lines (E\ < 40 K). Lines 
higher than P(5) (E\ > 40 K) are used to obtain the temperature 
of the warm component. Model 2 at t/t acc = 0.78 is used since it 
clearly shows the break between the two temperatures. 

How do the two temperatures evolve with time? As Table 2 
shows, the cold component is between 20 and 35 K with no sig- 
nificant evolution, except in Model 2, as the disk is building up. 
The warm component, however, increases with inclination and 
time. The derived warm temperature component ranges from 
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Fig. 12. Evolution of the 12 CO (red), 13 CO (blue) and C 18 
(black) absorption lines for Model 1 at i = 45°. The evolution 
of the absorption lines for Model 2 is shown at the bottom panel. 
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trum. All of the 13 CO lines are multiplied by 5 and C 18 lines 
by 25. 



< 100 up to ~ 520 K and traces the warming up of material 
in the inner region as the envelope dissipates. From the simula- 
tions, a > 400 K warm component is a signature of an evolved 
system. 

What is the origin of the two temperatures? This can be eas- 
ily shown by plotting the cold (<50 K) and warm (>100 K) H2 
column along the line of sight at various inclinations (Fig. 13). 
The particular model and time were chosen as an example; the 
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Table 2. NIR rotational temperatures derived from To > 10 3 lines for / = 45° and 75°. 
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Fig. 13. Integrated gas column density along the different view- 
ing angles for Model 2. The red and blues lines indicate the warm 
(T > 100 K) and cold gas (T < 50 K), respectively. The approxi- 
mate continuum optical depths, r c = 1 and 10, are also shown to 
indicate the amount material that is missed by properly solving 
the radiative transfer equation. 



general picture is the same independent of evolutionary model, 
but the warm component shifts to lower viewing angle (left) for 
earlier times. The figure indicates that the warm material can 
be probed through i > 25° viewing angles. The cold compo- 
nent probes the outer envelope where the gas temperature hardly 
changes as the system evolves. The warm component reflects the 
line of sight to the center through the warm inner envelope (the 
so-called 'hot core' ) or the disk surface where the gas is > 100 
K in all of our models. The cumulative gas column density in- 
dicates that the cold material is mostly absorbed at R > 100 AU 
while the warmer gas component is probing the inner few AU 



(i > 60°) up to -20 AU (i ~ 45°). The excitation temperature 
of the warm component reflects the density structure in the inner 
few AU. With the lower densities and warmer temperatures in 
Model 2, which decreases the effective continuum optical depth, 
higher warm temperatures can be observed as the absorption 
probes deeper to smaller radii than the other two models. 



5. Discussion 

We have presented the evolution of CO molecular lines based 
on the standard picture of a prestellar core collapsing to form 
a star and 2D circumstellar disk. The main aim of this work is 
to study the evolution of observables that are derived from the 
spectra such as rotational temperatures, line profiles, and veloc- 
ity fields. Specifically, signatures of disk formation during the 
embedded stage that can be obtained with the most commonly 
used molecule, CO, are investigated. 

5.1. Detecting disk signatures in the embedded phase 

How can one determine whether disks have formed in the ear- 
liest Stage 0? As discussed in Section 3.3, it is difficult to iso- 
late an embedded disk component in single-dish observations. 
The disk contribution to the observed flux within > 9" beams is 
only significant (>60%) for the 13 CO and C 18 J u > 9 lines at 
t/kee > 0.75. However, the absolute fluxes are too low to be de- 
tected. Thus, spatially resolved observations at sub-arcsec scale 
are needed. 

What are the signatures of an embedded rotationally sup- 
ported disk on subarcsecond scales? We have analyzed some 
of the simulated images during the Stage phase at a resolu- 
tion of 0.1 " as can readily be observed with ALMA. Figure 14 
shows the C 18 3-2 moment-one map during the Stage phase 
of Model 1 with a disk extending up to 20 AU at a resolution 
of 1, 0.3 and 0.1". Velocity gradients are present in all of the 
moment-one maps although it is weaker within al" beam. The 
right panel of the figure shows the spectroastrometry (Fig. 11) 
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Fig. 14. C 18 3-2 moment-one (left) and spectroastrometry 
(right) for Model 1 at the end of Stage convolved to a 1" 
(top), 0.3" (middle) and 0.1" (bottom) as shown by the solid 
circle. The Keplerian structure at this stage extends up to 20 AU 
(0.14" radius). The black dashed line in the right panel shows the 
v oc R~ 05 and the blue dashed line indicates the v oc R~ l relation. 

within the different resolutions. The velocity gradient within a 
0.5-1" beam exhibits a v oc R~ l relation since the inner envelope 
dominates the emission. Within the smaller beams, the velocity 
gradient is resolved into two components: v oc R~ l (blue dashed 
line) due to envelope and v oc R~ 05 (black dashed line) due to the 
stable disk. A v oc R~ 05 can be seen when the disk is marginally 
resolved, however, the derived stellar mass is far below the real 
stellar mass. It is very important that the disk is fully resolved in 
order to derive the stellar mass. With the C 18 lines, a resolved 
embedded disk in Stage would have the following features: 

- Elongated moment-zero map. 

- A transition between an infalling envelope to a rotating disk: 
a skewed moment-one map. 

- The Keplerian structure exhibits a clear v oc R~ 05 (Fig. 11 
and bottom of Fig. 14 ) velocity gradient perpendicular to 
the outflow. 

It is possible to search for rotationally supported disks in 
Stage sources with characteristics of Model 1 with the full 
ALMA array in less than 30 minutes for the C 18 2-1 transi- 
tion or up to 4 hours for the 6-5 transition. For Models 2 and 3, 
the disks have grown up to ~ 50 AU near the end of Stage and 
up to twice more massive, thus it will take less than 2 hours to 
observe them with ALMA in the 6-5 transition and less than 30 
minutes in the 2-1 and 3-2 transitions. Therefore, ALMA will 
allow us, for the first time, to test whether rotationally supported 
disks have grown beyond 10 AU by the end of the Stage phase. 

The question of when rotationally supported disks form 
is crucial to the physical and chemical evolution of accretion 
disks. In our evolutionary models, the disks have radii up to 
50 AU at the end of Stage phase and can be as small as 
20 AU. Dapp&Basu (2010) and Dapp et al. (2012) numeri- 
cally showed that disks do not grow beyond 10 AU during the 
Stage phase in the presence of magnetic fields due to the mag- 
netic breaking problem (see also Li et al. 201 1 ; Joos et al. 2012). 



This is consistent with the lack of observed rotationally sup- 
ported disks toward Stage objects so far (Brinch et al. 2009; 
Maury et al. 2010). On the other hand, there is a growing body 
of observational evidence for rotating disks in Stage I YSOs 
(Brinch et al. 2007; Lommen et al. 2008; J0rgensen et al. 2009; 
Takakuwa et al. 2012). The size of the rotationally supported 
disk depends on the initial parameters of the collapsing enve- 
lope. In our models, the sizes of Models 2 and 3 are consistent 
with the observed Keplerian velocity structure in Stage I sources 
(J0rgensen et al. 2009). 

An additional caveat is that our models lead to a Keplerian 
disk with a flattened inner envelope on similar scales as the 
stable disk. Hydrodynamical simulations, on the other hand, 
show a stable disk embedded in a much larger rotating flat- 
tened structure and with more turbulent structure (Brinch et al. 
2008; Kratteretal. 2010). Furthermore, the 3D simulations by 
Harsono et al. (2011) suggest that 50% of the embedded disk 
is sub-Keplerian due to interaction with the envelope. The pre- 
dicted moment 1 maps and spectroastrometry with the current 
models should be performed for these numerical simulations for 
comparison with ALMA observations. 

5.2. Probing the temperature structure 

The derived rotational temperatures from the submm and FIR 
lines (J u = 1-10) within > 9" beams do not evolve with time, 
in contrast with the continuum SED, and largely trace the outer 
envelope (Section 3.2). Even though the disk structure is hotter 
with time, the emission that comes from within that region is 
much smaller than the beam. The mass weighted temperature of 
the system which contributes to the overall emission is constant 
with time. The rotational temperatures also cannot differentiate 
between the three different evolutionary models. In addition, the 
peak of the CO SLED is constant throughout the evolution at 
J u = 4 for the 13 CO and C 18 0. The higher observed excitation 
temperatures for 12 CO and 13 CO in a wide range of low-mass 
protostars from submm and FIR data point to other physical pro- 
cesses that are present in the system such as UV heating and C- 
shock components that affect the / > 5 lines (Visser et al. 2012; 
Yildizetal. 2012). 

The NIR absorption lines toward embedded YSOs are com- 
plementary to the submm/FIR rotational lines. Both techniques 
probe the bulk of the cold envelope, but the NIR lines also probe 
the warm component because the lines start absorbing at the 
r ~ 1 surface deep inside the inner envelope. Consequently, the 
cold component present in the NIR lines should be similar to 
the component observed in the submm. Indeed, the typical ro- 
tational temperatures of 20-30 K found in the NIR models for 
C 18 are comparable to the values of 30-40 K found from the 
submm simulations. 

The warm up process is accessible through NIR molecular 
absorption lines. In Section 4.2, it is found that the warm com- 
ponent of the Boltzmann diagram evolves with time. The derived 
warm temperatures have a large spread depending on inclination 
and density structure in the inner few AU. A higher inner enve- 
lope and disk density correspond to a ~ 100 K warm component, 
whereas a more diffuse inner envelope has a warmer temperature 
that can be up to 500-600 K. Thus, as long as the inclination is 
known, the temperature of the warm component may give us 
a clue on the density structure of the inner few AU where the 
> 100 K gas resides. 

The predicted cold and warm temperatures can be compared 
with those found in NIR data toward Stage I low-mass YSOs. 
For the cold component, a temperature of 15 K has been found 
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for L1489 from 13 CO data (Boogert et al. 2002), and 10-20 K 
for Reipurth 50, Oph IRS43 and Oph IRS 63 from 13 CO and 
C 18 data (Smith et al. 2009; Smith et al. 2010), consistent with 
our envelope models. For the warm component, values ranging 
from 100 to 250 K are found for these sources from the iso- 
topolog data. These warm rotational temperatures are on the low 
side of the model values in Table 2 for the later evolutionary 
stages, and may indicate either a compact envelope (high in- 
ner envelope densities) and/or high inclination (which effectively 
limits the absorption lines to probe the inner envelope compo- 
nent). Independent information on the inclination is needed to 
further test the physical structure of those sources. 

6. Summary and conclusions 

Two-dimensional self-consistent envelope collapse and disk for- 
mation models have been used to study the evolution of molec- 
ular lines during the embedded phase of star formation. The 
non-LTE molecular excitation has been computed with an es- 
cape probability method and spectral cubes have been simu- 
lated for both the FIR and NIR regimes. The gas temperature 
is taken to be well coupled to the dust temperature as obtained 
through a continuum full radiative transfer code. With the avail- 
ability of spectrally resolved data from single dish submillimeter 
telescopes, the Herschel Space Observatory, VLT CRIRES and 
ALMA, it is now possible to compare the predicted collapse dy- 
namics with observations. We have focused the analysis on the 
13 CO and the C 18 lines since 12 CO lines are dominated by out- 
flows and UV heating. The main conclusions are as follows: 

- Spectrally resolved molecular lines are important in compar- 
ing theoretical models of star formation with observations 
and to distinguish the different physical components. The 
collapsing rotating envelope can readily be studied through 
the C 18 lines. Their FWHM probes the collapse dynamics, 
especially for the higher-/ (J u > 6) lines where up to 50% of 
the line broadening can be due to infall. 

- The derived rotational temperatures and SLED from submm 
and FIR pure rotational lines are found to be independent of 
evolution and do not probe the warm up process, in contrast 
to the continuum SED. The predicted rotational temperatures 
are consistent with observations for C 18 for a large sample 
of low-mass protostars. 

- The predicted 12 CO and 13 CO rotational temperatures and 
high-/ fluxes are lower than those found in observations of 
a wide variety of low-mass sources (Goicoechea et al. 2012; 
Yildiz et al. 2012, 2013). This indicates the presence of ad- 
ditional physical processes that heat the gas such as shocks 
and UV heating of the cavity walls (Visser et al. 2012). 

- The NIR absorption lines are complementary to the 
FIR/submm lines since they probe both the cold outer en- 
velope and the warm up process. Values obtained for the 
cold component are consistent with the observational data 
and models. For the warm component, the observed values 
are generally on the low side compared with the model re- 
sults. The high-/ NIR lines are strongly affected by radiative 
transfer effects, which depend on the physical structure of the 
inner few AU and thus form a unique probe of that region. 

- The simulations indicate that an embedded disk in both Stage 
(M env > Af*) and I (M env < M+ but M disk < M env ) does not 
contribute significantly (< 50%) to the emergent J u < 8 lines 
within > 9" beams. Higher-/ isotopolog lines have a higher 
contribution but are generally too weak to observe. The disk 
contribution is significantly higher within 1" beam and the 



evolution within such a beam indicate rotationally supported 
disks should be detectable in Stage I phase consistent with 
observations (J0rgensen et al. 2009). 

- Embedded disks during the Stage I phase are generally large 
enough to be detected with current interferometric instru- 
ments with 1" resolution, consistent with observations. On 
the other hand, high signal-to-noise ALMA data at ~ 0.1" 
resolution are needed in order to find signatures of embedded 
disks during the Stage phase. A careful analysis is needed 
to disentangle the disk from the envelope. 

- We have shown that the rotationally dominated disk can 
be disentangled from the collapsing rotating envelope with 
high signal-to-noise and high spectral resolution interfer- 
ometric observations (Section 3.4). The spectroastrometry 
with ALMA by plotting the velocity as function of peak po- 
sitions can reveal the size of the Keplerian disk. The 13 CO 
lines within interferometric observations may still be con- 
taminated by the infalling rotating envelope. Thus, more op- 
tically thin tracers are required. 

The three different collapse models studied here differ 
mostly in their physical structure and velocity fields on 10-500 
AU scales and illustrate the range of values that are likely to be 
encountered in observational studies of embedded YSOs. It is 
clear that the combination of spatially and spectrally resolved 
molecular line observations by ALMA and at NIR are crucial 
in determining the dynamical processes in the innermost regions 
during the early stages of star formation. A comparison between 
the standard picture presented here or hydrodynamical simula- 
tions and molecular line observations of Stage YSOs will fur- 
ther test the theoretical picture of star and planet formation. 
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Appendix A: Two dimensional RT grid 

The grid needs to resolve the steep density and velocity gradi- 
ents at the boundaries between adjacent components as shown 
in Fig. A.l. Improper gridding can lead to order-of-magnitude 
difference in the high-/ (/ > 6) line fluxes. The escape proba- 
bility code begins with a set of regularly spaced cells on a log- 
arithmic grid to resolve both small and large scales. Any cells 
where the abundance or density at the corners differs by more 
than a factor of 5, or where the temperature at the corners differs 
by more than a factor of 1.5, are split into smaller cells until the 
conditions across each cell are roughly constant. The full non- 
LTE excitation calculation takes about five minutes for a typical 
number of cells of 15 000 (Stage I) - 25 000 (Stage 0). 

Figure A.2 presents an example of the gridding in our mod- 
els. Such gridding is most important for early time steps where 
the outflow cavity opening angle is small. The refining ensures 
that the non-LTE population calculation converges and high-/ 
emission which comes from the inner region can escape. 



1500 f^[ 

1000 

500 



-500 

-1000 

1500 
1000 



< 







< 



500 



-500 



-1000 



10 4 




-1500 



-500 500 
R [AU] 



1500 



Fig. A.l. Gas density and velocity field for Model 3 at t/t^c - 
0.13 and 0.28 within 1500 AU. The density contours start 
from log(ft/cm -3 ) = 5.5 and increase by steps of 0.5 up to 
log(ft/cm -3 ) = 9. There are vertical and radial motions in the 
disk, which are not captured in this figure. 



Appendix B: FIR and submm lines 

Figure B.l shows the mass evolution of the evolutionary mod- 
els which indicates the times at which the models enter various 
evolutionary stages. Figure B.2 presents the C 18 9-8 lines for 
the three different evolutionary models convolved to a 9" beam. 
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Fig. A.2. Example of the gridding used in the Stage phase to 
resolve the small extent of the outflow cavity and to resolve the 
density and velocity gradient from the disk midplane to the enve- 
lope. The color contours are the same as in Fig. A.l. The struc- 
ture shown is for Model 1 at f/f acc = 0.13. 
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Fig. B.l. Evolution of the envelope, disk and stellar mass (top), 
effective stellar temperature (middle) and stellar luminosity (bot- 
tom) for the three different models as function of time (in units 
of t acc ). The solid circles show the time steps (roughly around 
t/kcc ~ 10" 3 , 0.1, 0.5, 0.75, 0.89, 0.99) used for rendering the 
molecular lines. The adopted time steps vary per model in or- 
der to cover properly the time when the model enters Stage 1 
(M env < M+) and when the M& = M env as indicated by the verti- 
cal dotted lines. 



Most of the line is emitted from the inner envelope region. The 
line is generally broader in later stages due to a combination 
of thermal, turbulent width and velocity structure. However, the 
width of the line depends on whether the emitting region is infall 
or rotation dominated. Figures B.3 and B.4 shows the J u =3, 6 
and 9 line profiles for the different models and how they change 
within different beams. The disk contribution to the observed 
lines within a 9" beam is presented in Fig. B.5 which is close to 
negligible for the low-/ lines. 
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Fig. B.4. As Figure3, but for Models 2 and 3. 
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AU as will discussed extensively in the next section. Figure C.2 
shows the inner 5 AU structure of the three different models. In 
both Models 1 and 3, the 10 4,5 cm" 3 density contour extends up 
to the line of sight of 45° due to the presence of compact disk and 
massive disk, respectively. On the other hand, the density struc- 
ture is flatter in Model 2. Thus, the NIR continuum emission is 
emitted through the high temperature region (mostly the red re- 
gion in the temperature plot) and encounters more material along 
the line of sight to an observer in Models 1 and 3. Consequently, 
the emission from the region close to the midplane of the disk 
will not contribute to the NIR continuum. Figure 13 presents an 
example of the radial contribution to the column density along 
a line of sight from the star. At viewing angles < 45°, most of 



I?- u "> riisn n o r ^i i a + n// u c + u the warm material is located at 10-40 AU from the star along 

Fig.B.2. C i5 9-8 line profiles convolved to a 9 beam for the , ,. _ . , t ., _ , tl . ,. , _° 



three models at t/t acc = 0.13 (red), 0.50 (green) and 0.96 (black). 

Appendix C: NIR molecular lines 

C.1. Inner disk 

Section 4 discusses the RADLite results with respect to the pre- 

The formation of NIR continuum and the region where absorp- dieted observables. In this section, we present the comparison 

tion happens depend on the physical structure in the inner few between RADLite and a simple calculation using the integrated 



the line of sight while most of the cold material is located > 50 
AU. The inner disk and high density region (< 1 AU) is only 
accessible through / > 75° in our models. 

C.2. Radiative transfer effects 
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Fig. B.3. Normalized C 18 and 13 CO line profiles as functions 
of evolution for J a = 3, 6 and 9 at f/f acc = 0.13 (red), 0.50 
(green) and 0.96 (black) for i = 5° orientation for Model 1. The 
lines are convolved to beams of 20" (top), 9" (middle) and 1" 
(bottom). 
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Fig. B.5. Disk contribution to total emission as a function of J u 
within a 9" beam (1260 AU at 140 pc) for the three different 
models. The different colors correspond to different time steps: 
t/t acc ~ 0.50 (red), 0.75 (blue), 0.96 (black). 



column density (column method, equation 3) to examine the ra- 
diative transfer effects in particular where the 4.7 yum continuum 
is formed and its affect on the simulated absorption lines. 

Figure C.4 illustrates the difference in line center optical 
depth between the formal solution of the radiative transfer equa- 
tion and the column density approach at i = 45° (Eq. 3). This 
particular model and time was chosen as an example that clearly 
shows the difference. The difference in To between the methods 
generally increases with J\ and are largest for J\ ~ 10. Figure C.3 
shows the 4.7 yum continuum flux due to dust only (i.e., no stellar 
contribution) as function of radius at various viewing angles. The 
continuum is point-like for low inclinations but it is generally 
off-center for high inclinations. Thus, larger tq differences be- 
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Fig. C.l. Example of the two-temperature fitting of the NIR 
lines. The column densities are derived from a curve-of- growth 
analysis performed on the simulated spectra using RADLite. The 
different symbols correspond to 12 CO (circle), 13 CO (triangles) 
and C 18 (squares). The solid lines are the combination of the 
cold and warm components. 
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Fig. C.2. Density structure in the inner 10 AU (starting from 
n(H2) = 10 4 ' 5 cm -3 ) and temperature structure (starting from 60 
K) of the three models. The three viewing angles are shown by 
dotted lines. 



tween the two methods are expected for larger inclinations since 
the absorption does not take place directly along the line of sight 
to the star. The implication is that the column method overesti- 
mates To of the high-/ absorptions as this method takes the full 
line of sight through the warm material into account. 

What are the implications for translating observed optical 
depths to column densities and excitation temperatures? The dif- 
ferences are largest for / > 5 lines which result in lower derived 
warm temperature and column densities in the full treatment. 
Consequently, the direct transformation from the observed opti- 
cal depth to column density requires additional information on 
where the absorption arises, i.e., the NIR continuum image. One 
way is to construct the physical structure and assess the NIR 
continuum image. Observationally, a comparison between inter- 
ferometric submm and NIR continuum should result in different 
continuum position if the NIR continuum arises from scattered 
light. Observationally driven constraint is likely more useful in 
practice since modelling of the NIR continuum requires a so- 
phisticated inner disk and envelope physical structure. For sys- 
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Fig. C.3. Normalized continuum flux without the stellar photo- 
sphere as function of distance from the center along one direc- 
tion. The figure is constructed from a NIR image of Model 1 
at t/tzcc = 0.7 as an example. The general trend is similar for 
^Aacc > 0.5 independent of evolutionary models. The normal- 
ized continuum is similar for the various evolutionary models 
for t/ta CC > 0.5. The continuum is shown to arise in the inner 
few AU region and to fall off very quickly for i < 60°. However, 
significant contributions from larger radii are expected for rela- 
tively high inclination. 
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Fig. C.4. Comparison of the line center optical depth obtained 
from a simple integration of the column density along the line of 
sight (circles) and the radiative transfer calculation (triangles). 
The P and R branches are shown in black and red, respectively. 
The horizontal line indicates the current observation limit with 
CRIRES at r ~ 3 x 10 -3 . The 13 CO lines are shown for Model 1 
aU/J acc = 0.76 for/ = 45°. 



terns with an off-centered NIR continuum, the derived column 
densities and temperature of the warm component are lower lim- 
its. 
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Fig. C.5. Comparison between LTE (squares) and non-LTE (tri- 
angles) line center optical depth obtained from RADLite for 
Model 3 at an inclination of 45°. The different colors correspond 
to the P (black) and R (red) branches. 
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Fig. C.6. Similar to Fig. C.5 for 12 CO (diamonds and circles) and 
13 CO (squares and triangles). 



the full non-LTE escape probability method (§3) or assuming 
LTE populations. The LTE level population is calculated us- 
ing the partition functions provided by the HITRAN database 
(Rothman et al. 2005). 

Figure C.5 compares LTE and non-LTE line center opacities 
of C 18 for Model 3 at an inclination of 45° (other isotopologs 
are shown in Figure C.6). Model 3 was chosen as an example; 
in general, the non-LTE effects are most apparent for the higher 
levels independent of isotopologs and evolutionary models. This 
reflects the findings of §3: lower pure rotational levels are more 
easy to thermalize due to lower critical densities. The difference 
between LTE and non-LTE is independent of inclination while 
radiative transfer effects do. Thus, it is more important to treat 
the radiative transfer correctly in the NIR to derive observables. 



C.3. Non-LTE effects 

The effects of non-LTE excitation in the vibrational ground state 
are studied by either using the level populations calculated with 



